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C/3 ■ Abstract 
O 

The Bootstrap method apphcation in simulation supposes that 
value of random variables are not generated during the simulation pro- 
^ I cess but extracted from available sample populations. In the case of 

. Hierarchical Bootstrap the function of interest is calculated recurrently 

CO \ using the calculation tree. In the present paper we consider the op- 

timization of sample sizes in each vertex of the calculation tree. The 
dynamic programming method is used for this aim. Proposed method 
CO . allows to decrease a variance of system characteristic estimators. 

CO I Keywords: Bootstrap method, Simulation, Hierarchical Calculations, 

Variance Reduction 

>< : 1 INTRODUCTION 



Main problem of the mathematical statistics and simulation is connected 
with insufficiency of primary statistical data. In this situation, the Bootstrap 
method can be used successfully (Efron, Tibshirani 1993, Davison, Hinkley 
1997). If a dependence between characteristics of interest and input data 
is very composite and is described by numerical algorithm then usually it 
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applies a simulation. By this the probabilistic distributions of input data 
are not estimated because the given primary data has small size and such 
estimation gives a bias and big variance. The Bootstrap method supposes 
that random variables are not generated by a random number generator 
during simulation in accordance with the estimated distributions but ones 
are extracted from given primary data at random. Various problems of this 
approach were considered in previous papers (Andronov et al. 1995, 1998). 

We will consider the known function of m independent continuos ran- 
dom variables Xi, X2, . . . , X^n : ^(Xi, X2, . . . , X^). It is assumed that distri- 
butions of random variables {Xj} are unknown, but the sample population 
Hi = {Xji,Xj2, . . . ,Xj„.} is available for each Xi,i = l,m. Here is the 
size of the sample Hi. The problem consists in estimation of the mathemat- 
ical expectation 



The Bootstrap method use supposes an organization of some realiza- 
tions of the values ^(Xi, X2, . . . , X^). In each realization the values of 
arguments are extracted randomly from the corresponding sample popula- 
tions Hi,H2, . . . ,Hm- Let be a number of elements which were ex- 
tracted from the population Hi in the l-th realization. We denote X(/) = 
(Xx X2j-2(/), • • • , Xjn^j„^{i)) and name it the l-th subsample. The estima- 
tor B* of the mathematical expectation Q is equal to an average value for 
all r realizations: 



Our main aim is to calculate and to minimize the variance of this esti- 
mator. The variance will depend upon two factors: 1) a calculation method 
of the function 0; 2) a formation mode of subsamples X(i). 

The next two Sections will be dedicated to these questions. In Section 4 
we will show how to decrease variance D Q* using the dynamic programming 
method. 



G = i?(/)(Xi,X2,...,X^). 



(1) 




(2) 
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2 THE HIERARCHICAL BOOTSTRAP METHOD 



We suppose that the function cp is calculated by a calculation tree. A root of 
this tree corresponds to the computed function cf) = <j)}^. It is the vertex num- 
ber k. The vertices numbers 1,2, ... ,m correspond to the input variables 
Xi,X2,... ,Xjn- The rest vertices are intermediate ones. They correspond 
to intermediate functions (f>m+ij 4'm+2, ■ ■ ■ , 4>k-i (see Fig.l). 




H.-, H 



Figure 1: Calculation tree 

Only one arc y„ comes out from each vertex u, u = 1, 2, . . . , fc. It cor- 
responds to a value of the function We suppose that = X^^v = 

1,2, . . . = = <Afe- 

We denote 1^ as a set of vertices from which arcs come into the vertex f , 
and Jy as a corresponding set of variables (arcs): i e ly <^ yi & Jy. It is clear 
that ly = (Z) ioi V = 1,2, ... ,m; yy = fv{Jv)j v = m + l,m + 2, . . . ,k. We 
suppose that a numbering of the vertices is correct: ii v <w than w ^ ly. 

Now, function value can be calculated by the sweep method. At the be- 
ginning, we extract separate elements Xij^(i), X2j2(i), • • • , from pop- 
ulations Hi,H2, . . . , Hjn, then calculate the function values (pm+i, 4'm+2, ■ ■ ■, 
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(j)k = (p successively. After r such runs the estimator Q* is computed ac- 
cording to the formula An analysis of this method was developed in the 
previous papers of authors. 

The Hierarchical Bootstrap method is based on the wave algorithm. Here 
all values of the function (j)^ for each vertex v = m + \,m + 2^...,k should 
be calculated all at once. They form the set Hy = {Yui, y^2j • • • > ^m^}, 
where is number of realizations (sample size) . Getting one realization Y^i 
consists of choosing value from each corresponding population Hi^i € /(w), 
and calculation of (p^ value. By this we suppose that a random sample with 
replacement is used when each element from Hi is choosen independendly 
with the probability l/n^. Further on, this procedure is repeated for next 
vertex. Finally we get l^i, ^^2, . . . , i^fon^ as values of the population H^- 
Their average gives the estimator Q* by analogy with formula ([2]). 

3 EXPRESSIONS FOR VARIANCE 

The aim of this section is to show how to calculate variance D Q* of the 
estimator ([2]). It is easy to see that in the case of Hierarchical Bootstrap 
the variance D Q* is function of sample sizes ni, n2, . . . , n^. 

In the previous papers of authors the variance D B* was calculated using 
the a;-pairs notion (Andronov et. al, 1996, 1998). Then, it was considered as 
continuos function of variables ni,n2, ■ ■ ■ ,nk, and reduced gradient method 
was used. But now we need other approach for the calculation of D Q*. 

We use Taylor decomposition of function 0(x) in respect to mean /i = 
{E Xi, E X2, . . . , E Xm)'- 

(/>(x) = </>(;x) + v^(/'(/i)(x-/x) + ^(x-/i)^v'0(/«)(x-/i)+O(||x-^I|3), (3) 

where \/'^(p{x) is the matrix of second derivatives of the function (p{x), 
||v|| is Euclidean norm of vector v. 

It gives the following decomposition: 
0(x)(/.(x') = (P\fi) + (/>(/i) 0(/i)(x - /i) + (/>(/x) 'A(/^)(x' - /x) + 
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+i<A(/x)(x'-^f v'<^(/^)(x'-^)+ (4) 

+ <P{^J){^ - /^)(x' - fif V + 

+0(||x'-/i[|3 + ||x-^||3). 

If X = (Xi,X2, . . . , is a random vector with mutual independent 
components {Xi}, E X. = fi, D Xi = af , then 

1 (9^ 

<A(X) = + - 5] ^.'^<A(/x) + ^(0(||X - (5) 

i=l * 

(i? = + '^(^) E ^'a^'^(^) + • • • • (6) 

z=l ^ 

Now we suppose that Xi and X^' are some values from sample population 
Hi = {Yii,Yi2, . . . ,Yin.}. Let Cj denote the covariance of two elements Yn 
and Yiii with different numbers / and 

Ci = Cov{Yii,Ya:\l^l'). (7) 

Let Xi and X'^ correspond to the elements Yn and Yn/ accordingly. Be- 
cause we extract Xi and X- from Hi at random and with replacement, then 
the event {/ ^ I'} occurs with the probability 1/nj. Then 



Cov{X„X'i 



Therefore 



Cov{Yii,Yii) = DYi with the probability 1/nj, 
CoviYii^Yiii) = Ci with the probability 1 — 1/nj. 

(8) 



Cov{Xi,Xi) = -DYi+(l--] a. (9) 
rii \ Ui,' 



If the values (j){X) and correspond to subfunction y„ = and 

the sample population H^, then formulas ([5D, (0) and ^ give 

^- = E f^'A.(/^.))'c'ot;(Xj,X;) + .... (10) 



ieiv 

Now we can get from 



Cov{MX),MX')) = —ol + f 1 - — ) a, (11) 



where the variance a'^ = D can be determined from (jlOp by Xi = X[: 

d 



E 

i£lv 



dx. 



(12) 



Finally we have 



Cov{Mx),Mx')) = E ( ^Mi^v 



I-—] Cov{X„Xi) 
rivJ 



+. 



(13) 



or 



Cov{Mx),Mx')) = J2 ( ^Mf^v'^ 



-+{l-l-]l]at + ... 

nr. 



n,, J n 



1 



71, 



+ 1- — i--a 
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(14) 



By this we suppose that variances erf , o"! i • • • i '^m of input random vari- 
ables Xi, X2, . . . , Xm are known. For example, it is possible to use estima- 
tors of these variances, calculated on given sample populations Hi, H2, ■ ■ ■ , Hm- 

If the vertex i belongs to the initial level of the calculation tree (it means 
that i = 1, 2, . . . , m) then (j)i{xi) = Xj, af is known value, Cj = 0. Therefore 

1 



Cov{Mx,i),Mx'J) = -^l 

m 



1,2, ... ,m. 



(15) 



Another covariances and variances are calculated recurrently in accor- 
dance with formulas (|10p . (|12p . (|14p . ()15p from vertices with less numbers to 
vertices with great numbers. Finally we get the variance D 0* of interest 
as the variance for root of the calculation tree: 

1 



DQ* = -{ai + (r - l)Cfc) + 



(16) 



where r = rik- 

As it was just mentioned, we will consider the variance D 0* as a function 
of sample sizes ni, 71-2, • • • , and denote it D{ni,n2, . . . , n^). Our aim is 
to minimize this function in respect to variables ni,n2, ■ ■ ■ ,nk by linear 
restriction, or, by other words, to solve the following optimization problem: 



minimize D(ni,n2, . . . , n^) 



(17) 
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by restriction 

aini + a2n2 + . • . + akUk < b, (18) 

where b, {oj} and {rii} are integer non- negative numbers. 

Now we intent to apply the dynamic programming method (Minox 1989). 

4 THE DYNAMIC PROGRAMMING METHOD 

Let us solve the optimization problem ()17p . (|18p . Our function of interest 
(/)(x) is calculated and simulated recurrently, using the calculation tree (see 
Section 2). In accordance to the dynamic programming technique, we have 
"forward" and "backward" procedure. 

During "backward" procedure, we calculate recurrently so-called Bell- 
man function ^y(a,z), v = 1,2, ... ,k, z = l,2,...,b, < a < 1. Let 
us consider the subfunction </>„(•), that corresponds to the vertex v. This 
subfunction directly depends on variables i G ly, which correspond to 
incoming arcs for the vertex v. Additionally (py depends on variables {xj} 
and {ni} from which there exists path from leaves to the vertex v of our 
calculation tree. Let By denote corresponding set of variables {xi} and 

{y^}■ 

Now we need to denote some auxiliary functions. Let us introduce the 
following notation 

i;iia) = aaf + {l-a)Cov{Xi,Xl), i = 1,2, . . . ,k, < a < 1. (19) 

Then we are able to write in accordance with (I13p : 

d 

Now we have from ([19l), (O and (HI]) 



V'.(a) = E (^Mf^v)y iaaf + (1 - a) \—af + (1 - —)Cov{Xi,Xl)\ ] 

^gj yOXi / I, Ifly fly J J 

E i^Uk^v))^ \U + ^) a? + (1 - a)(l - -)Cov{X„X[ 



so 

Ma) = E (^'^^(Z^")) (« + • (21) 

Note that it follows from pip and ()19p that our variance of interest ()16p 

is 

De* = MO). (22) 

Values ipv{oi) depend on the sample sizes for all i G B^. We will mark 
this fact as ipv{oi) = '4}y{a;ni : i £ By). 

Now we are able to introduce above mentioned Bellman functions: 

^>^,(Q;, z) = min V'i)(q;; nj : z € -Bt,), (23) 

where minimization is realized with respect to non-negative integer variables 
rii that are satisfied the linear restriction 

difii < z. (24) 



It is clear that optimal value of variance -D*(-) for the problem (|17p . (|18p 
is equal to $fc(0, h). 

Bellman functions $t,(a,z) are calculated recurrently from t> = 1 to 
V = k for < a < 1 and z = 1,2,..., 6. Basic functional equation of 
dynamic programming has the following form: 

$„(a, z) = min ^ (^'^^(^'')) (" + "^T^' •^*) ^^^^ 

where minimaiztion is realized with respect to non-negative integer variables 
riy and {zi : i E ly} that satisfy the linear restriction 

The initial values of <I>^(-) are determined with the tree leaves by formulas 
(dSI), (HHD and ([23]): 

^y{a, z) = aa^ + {1 - a) j . (^l = (27) 

[Z/ Cly\ 



= al(a + {l-a) ] V v = 1,2, . . . ,m, 

V [z/av\J 

where [z/a^] - integer part of number z/a^. 

Thus the "backward" procedure is a recurrent calculation of Bellman 
functions $t,(a, z) for v = 1,2, ... ,k, < a < 1, z = 1, 2, . . . , 6 by using 
formulas (|27|) . (|25p . Finally we get the minimal variance 

D* Q* = ^k{0,b). (28) 

To calculate the optimal sample sizes n J , n2 , . . . , n^, we should apply 
"forward" procedure of dynamic programming technique. At first, we find 
nl and {z* : i £ Ik} hy solving the equation 

<^,{0,b)=mmY,(^Mf^k)y '^i(^,z*] (29) 



ieik 

where minimization is realized by condition 

akul + ^z* < b. (30) 

Let I~^{v) = {uj : V £ ctfe = ^/i^t- Then we recurrently determine 

by analogy the rest n* and {z* : i € Iv} for v = k — 1, k — 2, . . . , m + 1: 

$i,(a7-i(^),z^) = min 2^ — (/>„(^^) ] aj-i(^) H j , z* j 

(31) 

by condition 

Moreover we put 

a^, = a7-i(^) + (1 - a7-i(„))/n* (33) 

Finally the optimal sizes n* for i = 1,2, ... ,m are determined by the 
following way: 

n* = [z*/a,]. (34) 
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